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We study the effect of spin-orbit coupling on both the zero-temperature and non-zero temperature 
behavior of a two-dimensional (2D) Fermi gas. We include a generic combination of Rashba and 
Dresselhaus terms into the system Hamiltonian, which allows us to study both the experimentally 
relevant equal-Rashba-Dresselhaus (ERD) limit and the Rashba-only (RO) limit. At zero tempera¬ 
ture, we derive the phase diagram as a function of the two-body binding energy and Zeeman field. In 
the ERD case, this phase diagram reveals several topologically distinct uniform superfluid phases, 
classified according to the nodal structure of the quasiparticle excitation energies. Furthermore, 
we use a momentum dependent SU(2)-rotation to transform the system into a generalized helicity 
basis, revealing that spin-orbit coupling induces a triplet pairing component of the order parameter. 

At non-zero temperature, we study the Berezinskii-Kosterlitz-Thouless (BKT) phase transition by 
including phase fluctuations of the order parameter up to second order. We show that the superfluid 
density becomes anisotropic due to the presence of spin-orbit coupling (except in the RO case). This 
leads both to elliptic vortices and antivortices, and to anisotropic sound velocities. The latter prove 
to be sensitive to quantum phase transitions between topologically distinct phases. We show further 
that at a fixed non-zero Zeeman field, the BKT critical temperature is increased by the presence of 
ERD spin-orbit coupling. Subsequently, we demonstrate that the Clogston limit becomes infinite: 

Tbkt remains non-zero at all finite values of the Zeeman field. We conclude by extending the quan¬ 
tum phase transition lines to non-zero temperature, using the nodal structure of the quasiparticle 
spectrum, thus connecting the BKT critical temperature with the zero-temperature results. 

PACS numbers: 03.75.Ss, 05.30.Fk, 47.37.-l-q, 74.25.Uv, 75.30.Kz 


I. INTRODUCTION 


Spin-orbit coupling, the interaction of a particle’s spin 
with its motion, is an essential ingredient in many quan¬ 
tum mechanical phenomena. In atomic physics, this ef¬ 
fect arises from the interaction between the electron’s 
magnetic moment and the magnetic field generated by 
the electron’s orbital motion, giving rise to the fine- 
structure splitting. In condensed matter physics, spin- 
orbit coupling leads to intriguing phenomena such as 
topological insulators [l|, the quantum spin-Hall effect 
Si and Weyl fermions S|. However, in these cases, the 
strength of the spin-orbit coupling is intrinsic and, more¬ 
over, the complex structure of the materials used is not 
always known, making theoretical modeling an arduous 
task. 

By contrast, ultracold atomic gases offer a versa¬ 
tile system in which parameters such as the interaction 
strength, the spin-imbalance, the dimensionality and ge¬ 
ometry can be freely adjusted i, making them ideally 
suited for quantum simulation of many-body systems. 
However, because the atoms used in ultracold gases are 
neutral, creating artificial spin-orbit coupling required 
the exploration of new techniques. More specifically, the 
use of two-photon Raman transitions was suggested the¬ 
oretically [a-laand shortly thereafter implemented for 
bosonic gases Q. Subsequently, spin-orbit coup ling was 
created in systems of non-interacting fermions |10l . Ill| . 


Recently, the interacting spin-orbit-coupled Fermi gas 
near a Feshbach resonance has also been realized 121 
and the formation of Feshbach molecules was investi¬ 
gated theoretically [l^. The type of spin-orbit coupling 
achieved in these systems is that of equal Rashba [ij 
and Dresselhaus [1^ strength (ERD), which up till now 
is the only form realized experimentally. 

These seminal experiments have sparked a wide range 
of suggestions for new experimental set-ups. Proposals 
to create spin-orbit coupling without the use of Raman 
dressing (which suffers from heating problems) include 
rf dressing with an atom chip [l^ and using ladder-like 
optical lattices m. Furthermore, many proposals have 
emerged for the creation of Rashba-only spin-orbit cou¬ 
pling, including the creation of degenerate dark states 
using tripod laser coupling [3 [13 generalizing the 
Raman scheme used in the aforementioned experiments 
[23 Hll . For a more complete overview of the experimen¬ 
tal achievements in this rapidly developing field, we refer 
to the following excellent review papers [22 h^ . 

The first theoretical studies of spin-orbit coupling in ul¬ 
tracold gases focused on the three-dimensional (3D) case, 
with either Rashba-only (RO) coupling [2^ - 1^ or ERD 
coupling [ 23113 . Recently, the two-dimensional (2D) RO 
case has also received wide attention [sil-l^. as well as 
the 2D ERD case [23 , in part due to the ex perimental 
creation of a 2D interacting Fermi gas 
relation to topological superfluids 


and by its 
40(1 . However, in 


2D, at non-zero temperature, a phase transition from a 
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quasi-condensate to a non-superfluid paired phase arises 
due to the Berezinskii-Kosterlitz-Thouless (BKT) mech- 
anism, which involves the unbinding of vortex-antivortex 
pairs 1^ . To capture the physics of this phenomenon 
it is essential to go beyond the saddle-point (mean-field) 
approximation and include fluctuations of the phase of 
the order parameter into the description (dsl - l^ . We per¬ 
formed this calculation for the 2D case with generic spin- 
orbit coupling, which was reported in a recent Letter [46l| . 
In the current paper, we will discuss the full mathemat¬ 
ical details of the aforementioned calculation, as well as 
present additional results for the zero-temperature case 
and for the BKT critical temperature. 

The remainder of this paper is divided in two parts: 
saddle point and fluctuations. In Sec. m we develop and 
discuss the saddle-point case. We start in Sec. Ill Al bv 
introducing the system Hamiltonian and setting up the 
functional-integral formalism, which we use throughout 
the paper. A derivation of the saddle-point thermody¬ 
namic potential is shown in Sec. Ill B1 Subsequently, in 
Sec. nsi we make a momentum dependent transforma¬ 
tion to the generalized helicity basis, which shows the 
emergence of a triplet component of the order param¬ 
eter. In Sec. umi we define the topologically distinct 
uniform superfluid phases of the system, based on the 
nodal structure of the quasiparticle excitation spectra. 
Finally, in Sec. Ill El we calculate the zero-temperature 
phase diagram as a function of the two-body binding en¬ 
ergy and Zeeman field. 

In Sec. imi we include fluctuations of the phase of the 
order parameter around the saddle point. Our main goal 
is to study the Berezinskii-Kosterlitz-Thouless (BKT) 
transition temperature, for which these fluctuations play 
a crucial role. In Sec. IIII Al we start by introducing 
the phase into our formalism, followed by a derivation 
of the effective action using the functional-integral adia¬ 
batic approximation in Sec. IIII Bl In Sec. IIII Cl the re¬ 
sulting effective action is then expanded up to quadratic 
order in the phase, leading to a phase-fluctuation part 
of the action. We conclude our calculation by deriving 
an analytic expression for the fluctuation thermodynamic 
potential in Sec. nmi Furthermore, we discuss the ef¬ 
fect of spin-orbit coupling on the sound velocities (Sec. 
IIII El) and on the vortex-antivortex structure of the sys¬ 
tem (Sec. IIII FI) . We then continue to Sec. IIII Gl in which 
we study the influence of spin-orbit coupling on the BKT 
critical temperature. Finally, in Sec. IIII HI we relate the 
zero-temperature results to the BKT critical tempera¬ 
ture, by discussing the evolution of the quantum phase 
transition lines at non-zero temperature. In Sec. IlYl we 
draw conclusions. 


II. FUNCTIONAL-INTEGRAL DESCRIPTION 
AT THE SADDLE-POINT LEVEL 

In this section, we set up the functional-integral for¬ 
malism at the saddle-point level and we discuss the 


ground states of the system, as well as the zero- 
temperature phase diagram. 

A. Setting up the formalism 

In this work, we use a functional-integral approach to 
calculate thermodynamic properties of the system. More 
specifically, we write the partition function as a sum over 
Grassmann fields ■i/' and '0, weighted by the exponential 
of the action functional S 

^ ^0’r,r,s^0’r,T,s ®^P [ ‘S'('0r,r,s: 0r,T,s)] ■ (1) 

Here, r = (x, y) and r indicate position and imaginary 
time, respectively, while s = {tj 4 -} denotes the spin-state 
(spin-up and spin-down) of the spin-1/2 fermions. The 
action can be related to the Hamiltonian density 'H via 
a Legendre transformation 


*5'('0r,r,S5 '^v,T,s) 



( 2 ) 

Our aim is to study a two-dimensional (2D) Fermi gas, 
where the spin-orbit coupling is a generic combination of 
Rashba and Dresselhaus terms. The Hamiltonian density 
of this system can be divided into three parts: % ='Hq + 
TLs + 'Hi- Note that for the remainder of this paper we 
use the units h = 2m = fcp = 1 . 

The first part of % is 

”^0 ^ ^ 0r,r,s [( Ms) *^s,s' ^ 2 ^z,ss'] 0r,r,sq 

s,s' 

( 3 ) 

corresponding to the single-particle sector. We work 
in the grand canonical ensemble, hence the use of the 
Lagrange multiplier /ig, which is interpreted as a spin- 
dependent chemical potential, thus allowing for popula¬ 
tion imbalance. Furthermore, denotes a Zeeman field 
perpendicular to the (x-y)-plane. Experimentally, this 
field corresponds to the intensity D of the Raman tran¬ 
sition between different hyperfine states: hz = — 0 / 2 . 
Finally, (Ji represents the Pauli matrix. 

The second part of T-L is 

— 2 0r,r,s i^kxO'y^ss' 0r,r,sq (4) 

s,s' 

corresponding to the spin-orbit terms. Here, we have de¬ 
fined a = (fR, -I- vu)/‘2 and 7 = (tiR — fD)/2, with tiR and 
wd being the Rashba and Dresselhaus coupling strength, 
respectively. Moreover, ki = —i{d/dl) is the momentum 
operator along the ^direction. 

The third part "Hi describes the interaction between 
fermions. We consider s-wave scattering, thus only tak¬ 
ing into account interaction between fermions in different 
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spin states. For a general two-body potential V[r — r'), 
the interaction term can be written as 

"Hi = y (r - r')V^r',T,il/’r,r,t- (5) 

However, in this work, we restrict ourselves to short- 
range interactions, which can be described by the contact 
potential V{y — v') = g5{Y — r'). 

Having set up the functional-integral formalism, we are 
ready to discuss the saddle-point approximation. 

B. Calculating the saddle-point thermodynamic 
potential 

In this section, our goal is to derive the saddle-point 
thermodynamic potential from the partition function 
shown in Eq. O- The difficulty in calculating the latter 
expression analytically lies in the fourth-order interaction 
term. A frequently used method to circumvent this prob¬ 
lem is to use the Hubbard-Stratonovich transformation 


Fourier transforming the action and applying the 
saddle-point approximation, the action can be written 
as 

-S' = i ^ + 

k,Wn. 

k,a;„ s ^ 

where ujn = (2n -I- l)7r//3 is the fermionic Matsubara fre¬ 
quency, k is the fermionic wave vector and is the chem¬ 
ical potential in spin state s. However, for the remainder 
of this paper, we choose g,s = and treat only a system 
with initial identical populations. In Eq. dS]) the fermion 
fields are ordered using the spinor notation 

= ('^k,t V'-k '0-k i) . (9) 

In this basis, the division into a quasiparticle-quasihole 
part and a spin-up/spin-down part is visible. The Hamil¬ 
tonian density appearing in Eq. m is 


exp 





VArVAr 


/ ArAr - - - \ 

dr (- lAr - I i 

^ 

where we denote r = {r, t}. This transformation de¬ 
couples the fourth-order interaction term in Eq. © into 
second-order terms, at the cost of inserting an additional 
functional integral over complex fields Ar,,- and A^^r- 
However, these fields have a physical meaning: they are 
interpreted as the fermion pair fields. In Eq. , we 
use the Bogoliubov channel in the Hubbard-Stratonovich 
transformation. It is also possible to use a different 
channel by using fields which represent the total density 
(Hartree channel) or the population imbalance density 
(Eock channel). However, for the description of superflu¬ 
idity, the Bogoliubov channel is the most natural. 

At this point, the fermionic functional integrals in the 
partition function can be calculated analytically, since 
the action has been made quadratic in the fermionic 
fields. It is not possible, however, to calculate the bosonic 
functional integrals analytically and one has to resort to 
approximations. Several ‘levels’ of approximation can be 
considered, starting with the crudest one: the saddle- 
point approximation. In this case, we assume that the 
order parameter is constant in time and space or, equiv¬ 
alently, that only its zero momentum component con¬ 
tributes 


X exp 


= vW^(q)<5..„.o|A|. (7) 

Hence, we assume that fermions pair at opposite mo¬ 
menta, ignoring the possibility of non-uniform superfluid 
phases (47l - l5fll | . In Eq. © , the factor is added to 

give |A| dimensions of energy, where /3 denotes inverse 
temperature and is the area of the 2D system. 


^(k) 


/Ck-h, -hl(k) 0 |A| \ 

-/i±(k) a + -|A| 0 

0 -|A| -a + h, -h^{k) ■ 

V |A| 0 -hlik) -^k-hj 

( 10 ) 


The emergence of the second term in Eq. © stems from 
the fact that operators have to be Weyl-ordered before 
they can be mapped onto Grassmann variables. This 
leads to additional terms due to the anti-commuting na¬ 
ture of the fermion operators. 

In the matrix shown in Eq. OT, = k^ — /r is the 
single-particle energy relative to the chemical potential, 
and h^(k) = h^ik) + ihy(k) is the spin-orbit field with 
components hx:{k) = —2jky and hy{k) = 2akx- It is 
noteworthy to mention that the Hamiltonian density can 
be written in terms of the Pauli matrices as 


'H(k) = Tx® (^kO-Q - /izO-z) - |A|Ty (g) ay 

-I- 2jkyTo ®ax- 2akxTz g) cr^, (11) 

where the Pauli matrices ai and Ti are associated with 
the spin-part and the particle-hole part, respectively. Us¬ 
ing this notation, it can be shown that quasiparticle- 
quasihole symmetry is preserved because 

Tx ao'H{k)Tx 1^1 ao =-'H*{-k). (12) 

Now, the fermionic functional integral can be calcu¬ 
lated, leading to the partition function 


Z = exp ( i ^ Tr In 




— /3 (*Wn + k^ — g) + 


(13) 












results in the saddle-point thermodynamic potential 
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Finally, performing the Matsubara summation, and re¬ 
placing the interaction strength g in terms of the two- 
body binding energy i?B via the relation 


1 

9 


dk 1 
(27r)2 2 fc 2 -h Eb ’ 


(14) 


J 


flsp - 


dk 


1 


(,„{ 


2 -1-2 cosh 


/3e(+)^ 


(k) I -bln 1 2 -I- 2 cosh /3e^ ^(k) - 


2 fc 2 -b Eb 


(15) 


In Eq. ()15|) . the momentum-dependent functions 

eW(k) = ^JEl + hl + |h^(k)|2 ± 2^Elhl + ^2|hx(k)|2, 

(16) 

represent the quasiparticle energies. Note that by us¬ 
ing Eq. o, we deliberately choose Eb to represent 
the two-body binding energy in the absence of spin-orbit 
coupling. In this way, Eb can be treated as an indepen¬ 
dent system parameter, which makes it easier to identify 
the direct effects of spin-orbit coupling. For a detailed 
study of the effect of spin-orbit coupling on the bound 
state energies resulting from interaction between spin- 
1/2 fermions, we refer to [l3l| . 

Having discussed the saddle-point thermodynamic po¬ 
tential, we investigate next the generalized helicity basis. 


C. The generalized helicity basis 

To gain insight into the effects induced by the pres¬ 
ence of spin-orbit coupling, it is instructive to transform 
the system to the generalized helicity basis, using a mo¬ 
mentum dependent SU(2) transformation. In this ba¬ 
sis, the Hamiltonian density in the non-interacting limit 
['H(k)]|^l^o becomes diagonal. The transformation ma¬ 
trix is given by 


/ Uk Uk 0 0 \ 

-v^. Uk 0 0 1 

0 0 u-k I ’ 

Vo 0 -u_k u_k/ 


(17) 


with the components of the eigenvectors being equal to 


Mk 


\ 


i 1 +__ 

21 ^hi + \h^ikwr 




v\^ = -e 


A 


( 18 ) 


1 


1 - 


y/hl -b |hj.(k)|2 / 


Here, the phase tpk is defined by h±{k) = |/i_L(k)|e*'^*‘. 
Applying the generalized helicity basis transformation to 


the full Hamiltonian density induces new components of 
the order parameter 


U'<nik)U 


0 ^bfr(k) A44^(k) 

A;^(k) A*^(k) -e^(k) 0 

0 -C4i(k)y 

(19) 


On the diagonal, the energies of the two generalized 
helicity bands are given by 


e^(k) = a - ^hl + \h^{k)\\ 

^4(k) = fk + y/hl -b |hu(k)|2. 


( 20 ) 


Furthermore, the order parameter in the generalized he¬ 
licity basis is now a tensor with components 

^•frfr('*) = “^T(k), 

^M(k) = As(k), /21') 

A^^(k) = -A^^(k), ^ ^ 

A^^(k) = A;^(k). 


Here, we identify, respectively, the singlet and the triplet 
component of the order parameter 


As(k) 

AT(k) 




^hl + \hAW 
^l(k) 


'Jhl + |/i_L(k )|2 


|A|, 

|A|. 


( 22 ) 


However, these two components are not independent, as 
they satisfy 


|As(k)p + |AT(k)p = |Ap. (23) 

At this point, it is important to emphasize that our 
system only has one order parameter. A, which is a com¬ 
plex scalar, because there is only s-wave interaction in 
the original spin-basis. It is only in the generalized he¬ 
licity basis that the order parameter can be decomposed 
into a singlet and triplet component and that a spino- 
rial structure arises [30| . The importance of the spinorial 
structure was also discussed recently in the context of 
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a repulsive Fermi gas M- For non-zero spin-orbit cou¬ 
pling, this triplet component cannot be fully suppressed 
by a Zeeman field, as it involves pairing between particles 
of equal generalized helicity. Hence, irrespective of the 
magnitude of the Zeeman field, the order parameter will 
always contain a triplet component. 

We thus see that transforming to the generalized helic¬ 
ity basis has effectively changed the local isotropic s-wave 
interaction into a non-local anisotropic interaction. The 
triplet component of the order parameter is not inher¬ 
ently present in the system, instead it is induced by the 
spin-orbit coupling. For a related discussion on the sin¬ 
glet and triplet components of the condensate fraction 
we refer to [52| |. 

Now that we have presented a detailed discussion of 
the generalized helicity basis, we continue by studying 
the various uniform superfluid phases of the system. 


D. Topologically distinct superfluid phases 

The Hamiltonian density of the system has four eigen¬ 
values: two quasiparticle energies ep^^(k), shown in Eq. 
(fTBl) . and two quasihole energies = —ep'^^(k). The 

structure of these excitation spectra can be used to dis¬ 
tinguish different uniform superfluid (US) phases. Let 
us look more closely at the quasiparticle branches. The 
(-l-)-branch is always gapped, whereas the (-)-branch can 
have nodes depending on the system parameters. In the 
language of the generalized helicity basis, we can write 
the second quasiparticle energy as 


i-US-0 



i-US-O^US-2 




FIG. 1. Overview of the topologically distinct uniform su¬ 
perfluid (US) phases, categorized by the nodal structure of 
the lowest quasiparticle energy branch. In Figs. (a),(b),(c) 
and (d), the values of p and |A| are held fixed while hz is 
increased. In Figs, (e) and (f), the values of |A| and hz are 
held fixed while p is decreased. Note that these phases only 
occur in the ERD case. 


4'^ (k) = V[^s(k)-|heff(k)|]2 + |AT(k)P. (24) 
Here, we have introduced the energy associated with 
the singlet channel Us(k) = |As(k)p, as well 

as an effective magnetic field hefr(k) = {h±{'k),hz) 
which is a combination of the spin-orbit and the Zee- 
man fields. This effective held can also be written as half 
the energy difference of the generalized helicity bands 
|heff| = [Cf|.(k) — 5 fi-(k)]/2, while the single-particle en¬ 
ergy is equal to the average energy of the helicity bands 

= Ktr(k) + Cfi(k)]/2. 

The lowest quasiparticle energy branch e4^(k) has 
nodes whenever the following two conditions are satis¬ 
fied simultaneously: 1 ) the effective magnetic field heff(k) 
and the singlet energy £'s(k) a-re equal in magnitude, and 
2 ) the triplet component AT(k) of the order parameter 
is zero. In the ERD case, AT(k) = 0 leads to kx = 0, 
which together with Es(k) = |/ieff(k)| gives the relation 
{ky — fj,)^ + |Ap = h^, yielding the possibility of having 
nodes at non-zero momentum. By contrast, in the RO 
case, or any other combination of Rashba and Dressel- 
haus terms, no nodes are present at non-zero momentum. 

The different possible phases in the ERD case are 
shown in Fig. [1] Let us describe these phases in more 
detail. When the Zeeman held is smaller than the order 


parameter {h^ < |A|), two phases can be discerned: 1) If 
the chemical potential /r > 0 , the system acquires an in¬ 
direct gap at non-zero |A|, and 2) if /i < 0 a direct gap at 
ky = 0 occurs. These phases are labeled i-US-0 and d-US- 
0, respectively. When the Zeeman held becomes larger 
than the order parameter > |A|), the quasiparticle 
spectrum acquires a nodal structure, depending on the 
value of the chemical potential. If /x > yTifAUjAp, the 
spectrum has two pairs of nodes (US-2 phase). When 
ImI < - |A| 2 , one pair of nodes is removed from 

the spectrum at fc = 0 and only one pair remains (US-1 
phase). Finally, when /i < —\/h1 — |Ap the hnal pair 
of nodes also vanishes at fc = 0 and the system becomes 
directly gapped (d-US-0 phase). To summarize, the dif¬ 
ferent topological phases can be classified as follows 


hz < |A| 
hz > |A| —> 


J/i > 0 i-US-0 

i/x < 0 d-US-0 

[■ fi > y'hl - |A|2 US-2 

I |/x| < _ |A|2 US-1 

U <-Vhl - |A|2 d-US-0 


(25) 


There are several effects induced by the nodal struc¬ 
ture of the quasi-particles energies at low temperatures 
(T << Tbkt)- First and foremost there is a dramatic 
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change in the momentum distribution of the system. For 
example, in the i-US-0 phase the momentum distribution 
is a smooth function, whereas in the US-2 phase disconti¬ 
nuities develop. Furthermore, both the isothermal com¬ 
pressibility and the spin susceptibility are non-analytic 
at the phase boundaries between the different uniform 
superfluid phases. This provides clear thermodynamic 
signatures of the quantum phase transitions. For more 
details we refer to [s^. Furthermore, in this paper, the 
emergence of nodes in the order parameter, when viewed 
in the generalized helicity basis, leads to anisotropies in 
the superfluid density tensor and sound velocities, as de¬ 
scribed in sections III C and III E. 

Having identified the topological nature of the uniform 
superfluid phases, we continue by discussing the ground 
state phase diagram. 


E. Zero-temperature phase diagram 

To find out which of the aforementioned superfluid 
phases occur, we investigate the zero-temperature phase 
diagram as a function of the two-body binding energy 
Eb and Zeeman field h^- More specifically, for a given 
(i?B, hi;)-point, we minimize the saddle-point free energy 
Fs-p = Hgp -f /in with respect to the order parameter 
|A|, while simultaneously solving the number equation 
i9r2sp/i9/i = — n in order to determine the chemical poten¬ 
tial /i. The resulting values of |A| and /i then determine 
which phase the system reaches, according to Eq. (051) . 
In Eig. [21 the resulting (Eb, h;;)-phase diagram is shown 
for several values of the spin-orbit coupling strength. 

Eigure|21a) shows the case without spin-orbit coupling. 
In this case, only the standard gapped superfluid phase 
(US-0) occurs, with a crossover between an indirect gap 
(i-US-0) at low binding energy and a direct gap (d-US- 
0) at large binding energy. At each value of the binding 
energy a critical Zeeman field exists at which a first or¬ 
der transition occurs from the US-0 phase to the normal 
phase. The driving mechanism behind this transition 
is the energy splitting between spin-up and spin-down 
fermions caused by the Zeeman field, which suppresses 
spin-singlet pairing at opposite momenta. The stronger 
the two-body binding energy between fermions, the larger 
the Zeeman field must be to break up the fermion pairs. 
This first order phase transition is also visible in Fig. 
H a), where the order parameter |A| is shown as a func¬ 
tion of hz, for several values of the binding energy. This 
figure reveals that the order parameter jumps discontin- 
uously to zero at a critical Zeeman field which depends 
on the value of Eb- 

Figures HDb) and (c) show the case of ERD spin-orbit 
coupling with wr/uf = wd/wf = 0.5 and wr/wf = 
vb/vf = 0.8, respectively (with % = r'F/2 and uf the 
Eermi velocity). These figures demonstrate the existence 
of the topological uniform superfluid phases US-1 and 
US-2 at zero temperature. An important difference be¬ 
tween the case without spin-orbit coupling and the ERD 


(a) v^=0v^=0 



(b) v^=0.5v^=0.5 



(c) y= 0.8 y= 0.8 



FIG. 2. Zero-temperature phase diagram, as a function of 
the two-body binding energy Eb and a Zeeman field hz, for 
different values of the spin-orbit coupling strength. In these 
figures we consider only equal-Rashba-Dresselhaus spin-orbit 
coupling (v = vr = ud): (a) u = 0, (b) v/ve = 0.5 and (c) 
v/vf = 0.8 (uf = vf/2 with uf the Fermi velocity). 


case is that in the former the system always transitions 
into the normal phase at high Zeeman held, whereas in 
the latter the system transitions into the US-1 phase. 
The origin of this difference lies in the triplet component 
of the order parameter that is induced by the presence 
of spin-orbit coupling, as demonstrated in section III Cl 
This triplet component cannot be suppressed by a Zee- 
man held, as it involves pairing between particles of equal 
generalized helicity. Hence, irrespective of the magnitude 
of the Zeeman held, the order parameter will always con¬ 
tain a triplet component and as a result will only become 
zero in the limit > oo. This behavior is shown explic¬ 
itly in Fig. EKb): at low values of the Zeeman held, the 
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'^>IAI v„=0v„=0 v„=0.8v^=0.8 



FIG. 3. Order parameter as a function of the Zeeman field h^. 
In the case without spin-orbit coupling, there is a first order 
phase transition at a given critical Zeeman field. In contrast, 
when (ERD) spin-orbit coupling is present, |A| only goes to 
zero in the limit —>■ oo. 


order parameter is approximately constant, while at high 
values it becomes suppressed but stays non-zero. Figures 
[H and [3] both coincide perfectly with the results of . 

Figures ^b) and (c) both show a triple point, in which 
the i-US-0 phase, the US-2 phase and the US-I phase 
meet. For low binding energy, the system undergoes two 
phase transitions with increasing Zeeman field: i-US-0 —>■ 
US-2 —US-1. At higher binding energy, the system un¬ 
dergoes only one phase transition: US-0 —US-1. With 
increasing binding energy, the region of the US-2 phase 
shrinks because the lower bound increases. This occurs 
because the gap in the quasiparticle spectrum of the i- 
US-0 phase increases with increasing binding energy, and 
thus a higher value of the Zeeman held is needed in order 

I 


to bridge this gap. Note also that this gap only disap¬ 
pears at fca; = 0 in momentum space. For all other values 
of k the gap is topologically protected by the presence 
of spin-orbit coupling. Here, we conclude our discussion 
of the saddle-point case and move on to include phase 
huctuations around the saddle point. 

III. PHASE FLUCTUATIONS AROUND THE 
SADDLE POINT 

In this section, we discuss the effects of phase fluctua¬ 
tions and its impact on the finite temperature phase di¬ 
agram, sound velocities and vortex-antivortex structure. 

A. Introducing the phase 

The scope of this work is to study the Berezinskii- 
Kosterlitz-Thouless (BKT) transition, in which phase 
fluctuations of the order parameter play a fundamental 
role. To introduce the phase, the complex field of the 
order parameter can be re-written as 

Ar,, = |Ar,,|e*®-. (26) 

Furthermore, we use the gauge transformation —>■ 

•i/ij. to make explicit the dependence of the action 

on the phase. 


Inserting Eq. (l26)) into the partition function of Eq. o after applying the Hubbard-Stratonovich transformation 
defined in Eq. m yields 


Z = 


'DfpT,r,s'Dfp 


r,r,s 


D|Ar,T|7^0r,T exp 


dr J dr [/o(r,T)+/s(r,T)-b/i(r,T)] 


(27) 


where the different parts of the action (single-particle, spin-orbit coupling and interaction) are given by 


4(1-, t) = i’T,T,a ( ^ - V? - M ■ 


7 fif) 7 1 

- *[Vr(0r,.)] • Vr - ^V?(dr,.) + -^NriOr.r)? ) 


/s(r,r) = 2J2s'^r,T,s 


sa 


d i dOj. 


dx 2 dx 


7l(r,T) — '^r,r,'|'V^r,T,4.Ar^r F V^r,r,4.V^r,T,tA^^r 


d i dOr T 

^^r.T^^r.T 
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V'r.T 


(28) 


Here, we used the symbol s ambiguously: s = {tj-i} when 
used as an index and s = ±I when used as a number. 
Since phase fluctuations provide the dominant contribu¬ 
tion to the physics in 2D, we ignore from this point on 
the contribution of amplitude fluctuations by assuming 
that the amplitude of the order parameter is constant: 
|Ar,r| = |A|. This still leaves three functional integrals 
to be calculated for the partition function described in 
Eq. dm). 


B. The adiabatic approximation 

In principle, the fermionic functional integral in Eq. 
dm) can be calculated exactly because the action is 
quadratic in the fermionic fields. However, we need first 
to transform the action to reciprocal space in order to 
eliminate the space-time derivatives. At this point in the 
calculation, the functional-integral adiabatic approxima¬ 
tion is used m. Ill, . We assume that the phase field 
Or.T varies slowly in space and time compared to a simi- 






















lar variation of the fermionic fields ipr,T,s and tjjr,T,s- As 
a result, for a given configuration of the phase field, the 
conhguration of fermionic fields can be coarse-grained by 
averaging over the ‘fast’ degrees of freedom. 

Given this approximation, we can Fourier transform 
the partition function in Eq. (1221) and calculate the 
fermionic functional integrals analytically. This proce¬ 


dure leads to 

Z = J exp (0r,r)] , (29) 

where the phase-only effective action is given by 


-^eff 


-Tr I In 

2 I 
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+ /3 

k.CJn 



(30) 


In this expression, /3 is the inverse temperature, is 
the area of the 2D system, k denotes the fermionic wave 
vector and a;„ = (2n-|- l)7r//3 is the fermionic Matsubara 
frequency. The matrix in Eq. (I30p has dimensions 4x4 
and can be written in terms of the 2x2 matrices 


M± 


± T/iz - Ck -hl(k)T/ii \ 
-h^{k) T h*^ Titon Cl) 

(31) 


and ID)± = ±zo-y|A|. 

The kinetic terms in Eq. dSU have been divided 
into phase-independent and phase-dependent contribu¬ 
tions, where we defined The phase- 

independent terms are ^k = k^ — /r and hz- The phase- 
dependent terms are + |:[Vr(0r,T)]^ and 

= —Vr(0r.r) ’ k. The spin-flip terms also contain 
a phase-independent contribution corresponding to the 
spin-orbit coupling field h^(k) = —2'yky -I- 2iakx and a 

fi d9 dB 

phase-dependent contribution — ict . 

In Eq. (|30|) , the two final terms emerge by taking into ac¬ 
count the anti-commuting nature of fermionic operators. 
This is done by Weyl-ordering these operators in second 
quantized form, before mapping them onto Grassmann 
variables. 


C. Expansion of the effective action up to 
quadratic order in the phase 

The exact calculation of the partition function shown 
in Eq. (l27l) requires knowledge of the eigenvalues of the 
matrix described in Eq. O- These eigenvalues are the 
solution of a quartic equation, and hence are too cum¬ 
bersome to be used in any analytic solution for the effec¬ 
tive action. Instead, we treat the spatial and temporal 
derivatives of the phase-field as a small perturbation and 
expand the effective action up to second order in terms 
of these derivatives. 

Let us first introduce a shorter notation Ak,w„(6*,90) 
for the 4x4 matrix in Eq. (sni)- By adding and 
subtracting the phase-independent part of the effec¬ 
tive action, we re-write the trace-log of this matrix as 


Tr{ln[Ak,a;„(0,0) -I- Fk(d, 90)]}, where 

-h\ 0 0 

4 (k) 0 0 

0 el(k) (hi)* ’ 

0 hi ei(k)/ 

(32) 

is the fluctuations part, with ei(k) = — Ck- Using 

this notation, the trace-log can be written as 

Tr{ln[Ak.^„(O,O) + Fk(0,90)]} 

=Tr{ln[Ak.<.„ (0, 0)]} + Tr{ln[I + A'i (0, O)Fk(0, 90)]}, 

(33) 

where the first term on the right-hand side is the saddle- 
point contribution, and the second term is due to phase 
fluctuations. To preserve the readability of the paper 
we refer the details of the explicit calculation of ^eff to 
appendix Here, we just show the final expression for 
the effective action Ses = Ssp + Sr, where S'sp is the 
saddle-point contribution and 



(34) 


is the fluctuation action with v = {x,y}. 

Due to the presence of anisotropic spin-orbit coupling, 
the superfluid density tensor has unequal components 
Pxx ^ Pyyi except in the isotropic Rashba-only case (as 
well as the Dresselhaus-only case), where p^x = Pyy Us¬ 
ing a simple scale transformation, the action in Eq. (1341) 
can be written in a form that is equivalent to the action 
in the case without spin-orbit coupling: 

-Sfl = iy dry dr -k Ps [V(0r,r)]^ 

(35) 

where the effective superfluid density is now given by 
Ps = yjPxxPyy The exact expressions for the compress¬ 
ibility A and the superfluid density components PxxiPyy 
are given in appendix 



Fk(0,90) = 


/ 4(k) 

- KY 
0 

V 0 
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D. The fluctuation thermodynamic potential 


Finally, the fluctuation thermodynamic potential can 
be derived from the action given in Eq. (1351) . When 
calculating the functional integral defined in Eq. (0^1) . 
one has to be careful not to double count the fields, since 
is a real field, hence Here, q is 

the bosonic wave vector and vora = 21 : 11 X 1 /fi is the bosonic 
Matsubara frequency. A way to circumvent this problem 
is to write the partition function over half the total q- 
space: 


2 ,= n 

<?x>0 




de, 


* 


X exp 


/ 

- + 
\ qx>0 




Now the known result of a Gaussian bosonic functional 
integral can be applied, which leads to 


Zfl = exp 


q,rom / 


After performing the bosonic Matsubara summation, and 
using Zfl = exp(—we arrive at 


t 2 j-oo 

= <^<->"( 1 --'’"“)- (V) 

In 2D, the integral in this expression can be calculated 
analytically, yielding 


= -C(3) 


a 


27r/33 ’ 

which corresponds to the additional pressure 


(38) 


5P=^—T^, (39) 

2t: Ps 

created by exciting collective (sound) modes, to be dis¬ 
cussed in the next section. 


E. Sound velocities 

In the action shown in Eq. (l34l) . the presence of spin- 
orbit coupling results in an anisotropic superfluid density 
tensor (except in the isotropic RO case). As a result, the 
sound velocities of the system become anisotropic. In 
the case without s pin-or bit coupling, the sound velocity 
is equal to c = i/ps/A, where A is the compressibility 
and Ps is the superfluid density. In the present case, the 
sound velocities in the x- and y-direction are given by 



£'^=0.05 v/?=V£)=0.8 



FIG. 4. Sound velocities as a function of the Zeeman field 
hz, for a given value of the binding energy Ab/Af = 0.05 
and spin-orbit coupling strength vr/vf = vd/vf = 0.8. The 
sound velocities are shown at three different temperatures: 
T/Tf=0, T/Tf=0.02 and T/Tf=0.04. The presence of ERD 
spin-orbit coupling induces an anisotropy in the sound ve¬ 
locities. Moreover, they are sensitive to the quantum phase 
transition between the US-2 and US-1 phase [4^. Notice that 
the vertical axis does not start at zero. 


respectively. The anisotropic sound velocities of the 
system are shown in Fig. |4] as a function of the Zeeman 
field hz, at different temperatures. The binding energy 
is held fixed at Es/Ep = 0.05 and we consider ERD 
spin-orbit coupling with ur/uf = vu/vp = 0.8. This 
means that we follow a vertical line through the phase 
diagram in Eig. HKc). Along this line with increasing h^, 
the following uniform superfluid phases are encountered: 
i-US-0,US-2 and US-1. The transition points between 
these phases are indicated on the abscissa of Fig. |4l In 
this figure, we see that the sound velocities are sensitive 
to the presence of the quantum phase transition between 
the US-2 and US-1 phases: at this transition point, both 
sound velocities show a distinct cusp. On the other hand, 
the sound velocities do not show such behavior at the 
transition point between i-US-0 and US-2. 

The reason for this difference in sensitivity is that at 
the US-2-to-US-l transition, two nodal Dirac quasiparti¬ 
cles with opposite topological charges annihilate at zero 
momentum, i.e. in the long-wavelength limit. Because 
sound waves are low-energy and long-wavelength exci¬ 
tations, they tend to be sensitive to this transition at 
k = 0. The i-US-O-to-US-2 transition, however, can be 
understood (when approached from the US-2 side) as 
the annihilation of two Dirac quasiparticles with opposite 
topological charges at non-zero momentum. Therefore, 
the sound velocities are much less sensitive to this quan¬ 
tum phase transition. 

As expected, the sharpness of the cusps in the sound 
velocities tends to soften when temperature is increased. 
Note also that in the limit of —>■ 0, both sound ve¬ 
locities converge to the known limit Cx = Cy = vp/\/2 
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of the case without spin-orbit coupling. This shows that 
in the ERD case, the effect of spin-orbit coupling can be 
gauged away in the absence of the Zeeman field h^. Now 
that the sound wave excitations have been discussed, we 
proceed by analyzing the topological excitations. 


F. Vortex-antivortex structure 

The anisotropy of the superfluid density also has reper¬ 
cussions on the vortex and antivortex structure of the 
system. The vortex solutions are found by extremizing 
the fluctuation action (1351) in the static case (r = 0). This 
leads to Laplace’s equation V^(0r) = 0, which has sin¬ 
gular, vortex-like solutions that are given by 9±{x,y) = 
± arctan(y/a;). However, since the original fluctuation- 
action is given in Eq. (l34)) , we still need to transform the 
solution in order to reflect the rescaling that was per¬ 
formed in order to transform Eq. (IMl) into Eq. (1351) . 
This leads to the solution 


^±( 2 ;, y) = ± arctan 



(41) 


for a vortex(-|-)/antivortex(—) located at {x = Q,y = 
0). This shows that the vortices)V) and antivortices(A) 
present in the system become elliptical, rather than cir¬ 
cular, in the presence of ERD spin-orbit coupling. For 
the isotropic Rashba-only case, vortices remain circular. 

The general solution for a vortex-antivortex (VA) pair 
is found by taking the sum of two vortex-solutions with 
opposite winding number, leading to 


0vA(a:, y) = arctan (^ -2 ^ 2 ) ’ (^2) 

where x = x/p, y = yp and a = a/p, with p = 
[pxx/PyyY^'^- The parameter a indicates the distance 
between the cores of a vortex and an antivortex located 
at (x = —a,y = 0) and {x = a,y = 0), respectively. 
Plots of the solution given by Eq. (|42|) are shown in Fig. 
[SKa) for the RO case and in[SKb) for the ERD case. The 
parameters used in this figure are Eb/Ef = 0.01 and 
hz/Ep = 0.2. We chose these parameters to enhance vi¬ 
sualization, as the ratio of Pyy/pxx is larger for smaller 
binding energy. 

The emergence of elliptic vortices and the structure 
of the VA pairs in a 2D Fermi superfluid constitute im¬ 
portant signatures for the experimentally relevant ERD 
case. These signatures could be detected during a time- 
of-flight expansion of the trapped system, or via Bragg 
spectroscopy, which is also sensitive to the direction of 
rotation of the supercurrents (i^ . 

Having discussed the emergence of vortex-antivortex 
pairs, we investigate next the BKT transition tempera¬ 
ture, where vortex-antivortex unbinding occurs. 



FIG. 5. Vortex-antivortex structure for a 2D Fermi gas in (a) 
the Rashba-only case and (b) the equal Rashba-Dresselhaus 
case. In the presence of ERD spin-orbit coupling, the vortices 
and antivortices become elliptic. In the isotropic Rashba- 
only case, the vortex structure remains circular. The param¬ 
eters used are hz/Ep = 0.2, Eb/Ep = 0.01, T ~ Tbkt, with 
v-r/vf = 1 in (a) and vr/vf = vo/vf = 1 in (b). [i^ . 


G. Berezinskii-Kosterlitz-Thouless critical 
temperature 


In this section, we turn to the analysis of 
the Berezinskii-Kosterlitz-Thouless critical temperature 
(Tbkt)- In order to determine this transition tempera¬ 
ture, three equations need to be solved self-consistently: 
(1) the order parameter equation, determined by the 
condition i90sp/9|A| = 0, (2) the number equation 
—dflsp/dp = n, and (3) the generalized Kosterlitz- 
Thouless condition Tbkt = f Ps(7bkt), where ps = 
y^PxxPyy The first two equations define the saddle-point 
(mean-field) ‘transition’ temperature Tmf, at which the 
system changes from its normal phase with |A| = 0 to 
its paired phase with |A| yf 0. However, the transition 
to a true superfluid phase (quasi-condensate) occurs at 
Tbkt < Tmf, which is greatly affected by phase fluctua¬ 
tions. Determining this critical temperature requires the 
simultaneous solution of all three aforementioned equa¬ 
tions. 

Here we note that although the equation Tbkt = 
§Ps(Tbkt) was derived originally in the framework of 
the phase-only XY-model, it is still justified to use it 
in the present case with spin-orbit coupling. The rea¬ 
son being that the form of the phase-only action with 
spin-orbit coupling (IMl) . can be re-scaled to the phase- 
only action without spin-orbit coupling (1351) . All the 
effects of spin-orbit coupling are then contained in the 
effective superfluid density ps = y/PxxPyy- Because the 
resulting action is still of the same form as the action 
for the XY-model, the system with and without spin- 
orbit coupling fall within the same universality class, 
when the phase transition is continuous. As such, a sim¬ 
ple renormalization group analysis leads to the relation 
Tbkt = fPs(TBKT), where ps is the effective superfiuid 
density, defined above. 

In Fig. ini the BKT temperature is shown as a func¬ 
tion of the three main system parameters: the spin-orbit 
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coupling strength in (a) and (b), the binding energy in 
(c) and (d) and the Zeeman field in (e) and (f). Both the 
equal-Rashba-Dresselhaus (ERD) case and the Rashba- 
only (RO) case are shown, in the left and right column, 
respectively. The parameters vr, ud, Er, are the main 
energy scales of the system, and as such their relative 
magnitude will significantly affect Erkt, as is discussed 
in detail below. 

We begin by looking at the effect of the spin-orbit 
coupling strength on the critical temperature. In Fig. 

(a) and (b), Trkt is shown as a function of the spin- 
orbit coupling strength in the ERD case (ur = wd = 'c) 
and the RO case, respectively, for several values of the 
Zeeman field hz- Let us hrst consider the case without 
Zeeman field (blue filled circles). In this situation, the 
spin-orbit strength has no influence on the critical tem¬ 
perature in the ERD case. By contrast, in the RO case, 
the critical temperature is lowered with increasing spin- 
orbit coupling. The origin of this difference lies in the 
nature of the spin-orbit coupling itself. In both the ERD 
and the RO case, the spin-orbit Hamiltonian can be seen 
as introducing a gauge held, which leads to orbital mo¬ 
tion of the fermions. This makes pairing with opposite 
momenta harder (orbital frustration), thus suppressing 
superhuidity. However, in the ERD case, this gauge held 
can be removed by a gauge transformation, provided that 

hz = 0. 

Subsequently, we consider the case of non-zero Zee- 
man held. In the absence of spin-orbit coupling, it is 
known that Zeeman helds suppress singlet pairing and 
thus superhuidity. As shown in Fig.® (a) and (b), the 
introduction of spin-orbit coupling increases the critical 
temperature compared to the case wr = wd = 0. The 
reason for this is that a triplet pairing component is in¬ 
troduced by the presence of spin-orbit coupling, as was 
demonstrated in section Hi Cl where the generalized helic- 
ity basis was discussed. This triplet component cannot 
be suppressed by the Zeeman held. 

In the ERD case, Trkt rises monotonically and con¬ 
verges to its limiting value without spin-orbit coupling 
and Zeeman held in the limit n —^ oo. This occurs be¬ 
cause in the limit n —?► oo, hz becomes negligible and the 
ehect of ERD spin-orbit coupling can be gauged away. In 
the RO case, however, the situation is more complex. For 
hz = 0, the maximum lies at ur = 0. For low values of 
hz, Jbkt rises for small values of ur, reaches a maximum, 
and then decreases for increasing values of ur. For high 
values of hz , Ibkt is a monotonically increasing function 
of tiR that never exceeds the limiting value of Erkt at 
fR oo. 

In Fig. a) and (b) we have encircled two points. 
These points indicate that the nature of the critical tem¬ 
perature is different there. All other points indicate 
Tbkt, meaning that the generalized Kosterlitz-Thouless 
condition Trkt — is positive for T > Trkt and 

negative for T < Trkt- However, at these special lo¬ 
cations, the situation is different: Tbkt ~ is still 

negative for temperatures below this point, but at higher 


temperatures, the order parameter amplitude | A| is zero, 
meaning that the mean-held temperature Tmf has been 
reached. 

In Fig. [S] (c) and (d) the Berezinskii-Kosterlitz- 
Thouless critical temperature is shown as a function of 
the binding energy E-q , for different values of the Zeeman 
held hz, for the ERD and the RO cases, respectively. 
In both hgures, the cases with spin-orbit coupling are 
indicated by symbols, whereas the cases without spin- 
orbit coupling are indicated by dashed lines. Notice that 
Tbkt increases monotonically with increasing binding en¬ 
ergy, regardless of the value of the spin-orbit coupling 
strength. This is as expected, since a larger binding en¬ 
ergy leads to more strongly bound fermion pairs, thus 
strengthening superhuidity. In the asymptotic limit of 
inhnite binding energy, Tbkt converges to its limiting 
value: Tbkt/Tf = 0.125, regardless of the strength of the 
Zeeman held or spin-orbit coupling. This limiting value 
can be derived from the generalized Kosterlitz-Thouless 
condition: Tbkt = Tps/2. The superhuid density can at 
most reach half the total density, where the latter in 2D 
and in units h = 2m = fcp = 1 is given by 1 / 27 r. This then 
leads to Tbkt < 1/8. However, when the pairs become 
bosonic in nature with pair size ^pair much smaller than 
the interparticle spacing fcp logarithmic corrections due 
to boson-boson interactions may actually reduce Tbkt, 
asymptotically [ssLI^. 

We deliberately chose not to include these logarith¬ 
mic corrections to the BKT-critical temperature in the 
present treatment. Including these corrections involves 
a more complicated calculation, which lies beyond the 
scope of the present paper. It is important, however, to 
make sure that we do not reach the regime where log¬ 
arithmic corrections are relevant. To this end, we have 
calculated the pair size for the case without spin-orbit 
coupling (Eq. (19) in HI) and compared it with the 
average interparticle distance, given by kp^. In our pa¬ 
per, the latter quantity equals 1 because of our choice 
of units. In our case, without spin-orbit coupling, for 
hz/Ep = {0.01,0.21,0.41} and Ep/Ep varying between 
0.03 and 0.35, the pair size varies between 3.32 and 0.84. 
Thus, at these values of the binding energy, the regime 
where /'pair << kp^ has not been reached, validating our 
approach. 

In the ERD case, we see again that without a Zeeman 
field, the situations with and without spin-orbit coupling 
(blue filled circles and blue dashed line, respectively) are 
exactly equal. When hz becomes non-zero, however, the 
effect of spin-orbit coupling on the critical temperature 
is clear, as Tbkt is significantly increased. The effect is 
greater where hz is larger. Notice also that while Tbkt is 
a smooth curve in the presence of spin-orbit coupling, a 
jump occurs in the case without spin-orbit coupling, due 
to a first order phase transition. 

The RO case, shown in Fig. [IKd) is qualitatively sim¬ 
ilar to the ERD case. The main difference is that in the 
ERD case, at fixed hz, spin-orbit coupling increases Tbkt 
for all values of the binding energy compared to the case 
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FIG. 6. Berezinskii-Kosterlitz-Thouless critical temperature, for both the equal-Rashba-Dresselhaus (ERD) case and the 
Rashba-only (RO) case, shown in the left-hand and right-hand column, respectively. We have plotted Tbkt as a function of 
the three main system parameters: the spin-orbit coupling strength (a) and (b), the two-body binding energy (c) and (d) and 
the Zeeman field (e) and (f). All points shown indicate the BKT critical temperature, except for the two encircled points in 
(a) and (b), which indicates the mean-field critical temperature Tmf. 


without spin-orbit coupling. By contrast, in the RO case, 
this increase with respect to zero spin-orbit coupling oc¬ 
curs only at low binding energy, while at large binding 
energy Tbkt decreases, because of the aforementioned 
orbital frustration. The binding energy at which Tbkt 
crosses the critical temperature without spin-orbit cou¬ 
pling becomes larger with increasing Zeeman field. 

Finally, we study the critical temperature as a func¬ 
tion of the Zeeman field shown in Fig. HJe) and (f). 


The main point of these two figures is the fact that the 
Clogston limit (i.e. the critical Zeeman field at which 
the order parameter becomes zero) becomes infinite: 
the BKT critical temperature only becomes zero in the 
limit hz ^ oo. This finding is consistent with the zero- 
temperature behavior, as discussed in section ing There, 
it was shown that the US-1 phase survives all finite val¬ 
ues of the Zeeman held. Because our huctuation theory 
only describes the 2D system at non-zero temperature 
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it is comforting that the T —>■ 0 limit is recovered. In 
the following section, we connect the zero-temperature 
phase diagram and its quantum phase transition lines to 
the Berezinskii-Kosterlitz-Thouless critical temperature. 


H. Summarizing 3D phase diagram 

In Fig. [71 we show the finite temperature three- 
dimensional phase diagram as a function of the two-body 
binding energy and Zeeman field. This figure connects 
the zero-temperature phase diagram to the Berezinskii- 
Kosterlitz-Thouless critical temperature, in the ERD 
case. The zero-temperature phase diagram is plotted in 
the (FIb, h^)-plane. This phase diagram contains three 
uniform superfluid phases: i-US-0 (red circles), US-2 
(blue diamonds) and US-1 (green six-point stars). The 
BKT critical temperature as a function of Ub and is 
given by the red translucent surface. This surface only 
reaches T = 0 in the limit —>■ oo (for all values of E-q), 
as mentioned in the previous section. 

Finally, we have extended the quantum phase tran¬ 
sition lines to non-zero temperature. The ‘transition’ 
between i-US-0 and US-2 is given by the blue surface, 
and is defined as the finite temperature locus where the 
quasiparticle nodes of the US-2 phase merge at non-zero 
momentum, leading to the i-US-0 phase. We find that 
this ‘transition’ has a small temperature dependence: be¬ 
tween T = 0 and Tbkt this transition line is displaced 
by an amount on the order of O.lh^. Similarly, the ‘tran¬ 
sition’ between US-2 and US-1 is given by the green sur¬ 
face, which has a negligible dependence on temperature. 

To the best of our knowledge, no accurate theory 
exists to calculate the behavior of these lines at non¬ 
zero temperature. Here, we obtain these lines based 
on the quasiparticle description at finite temperatures. 
This can be regarded as an extrapolation between the 
zero-temperature results and the Berezinskii-Kosterlitz- 
Thouless temperature, both of which are known to be 
qualitatively accurate. Therefore, in the absence of a 
full theory, it is difficult to claim that these ‘transition’ 
lines are truly thermodynamic phase boundaries between 
topologically distinct phases, rather than crossover lines 
ending at topological phase transition lines at T = 0. 
We call them ‘transition’ lines because finite tempera¬ 
ture topological invariants can be defined on either side 
of these lines, provided that a quasiparticle component 
is still present. 


IV. SUMMARY AND CONCLUSIONS 

In this work, we have investigated the effects of spin- 
orbit coupling on both the zero-temperature and non¬ 
zero temperature behavior of a 2D Fermi gas with at¬ 
tractive interactions. We used a generic combination 
of Rashba and Dresselhaus terms, which allowed us to 
study both the equal-Rashba-Dresselhaus (ERD) and the 


Rashba-only (RO) limit. 

In the first part of the paper, we focused on results at 
the saddle-point level. Starting from the partition func¬ 
tion, we derived the thermodynamic potential within the 
saddle-point approximation. By minimizing this thermo¬ 
dynamic quantity, the zero-temperature phase diagram 
was obtained, as a function of the two-body binding en¬ 
ergy and Zeeman field. In the ERD case, we identi¬ 
fied several topologically distinct uniform superfluid (US) 
phases, classified according to the nodal structure of the 
quasiparticle energy bands. We distinguished between 
the uniform superfluid phase with zero, one or two pairs 
of nodes (US-0, US-1 and US-2 phases). We found that 
at any non-zero value of the spin-orbit coupling strength, 
the system is always in a US phase. More specifically, the 
US-1 phase survives at any large finite value of the Zee- 
man field. We identified this behavior by making a mo¬ 
mentum dependent transformation to the helicity basis, 
which diagonalizes the non-interacting Hamiltonian. In 
this basis, we showed that the order parameter acquires 
a triplet pairing component, which cannot be suppressed 
by a Zeeman field. 

In the second part of the paper, we focused on fluc¬ 
tuations around the saddle point. By expanding the 
action up to second order in the phase, the total ther¬ 
modynamic potential was written as a saddle-point and 
a fluctuation contribution. The latter contribution was 
rescaled to the corresponding action without spin-orbit 
coupling. We found that the superfluid density becomes 
anisotropic, due to the presence of spin-orbit coupling 
(except in the isotropic RO limit). We showed further 
that the anisotropic sound velocities are sensitive to the 
quantum phase transition between the US-2 and US-1 
phases, and that vortices and antivortices become el¬ 
liptical instead of circular. Subsequently, we studied 
the Berezinskii-Kosterlitz-Thouless critical temperature 
(Tbkt), by simultaneously minimizing the free energy, 
solving the number equation and satisfying the general¬ 
ized Kosterlitz-Thouless condition. Our three main find¬ 
ings were as follows. (1) Without the presence of a Zee- 
man field, ERD spin-orbit coupling can be removed by 
a gauge transformation, hence Tbkt remains unchanged 
compared to the case without spin-orbit coupling. In 
the RO case, however, increasing the Rashba coupling 
strength in the absence of a Zeeman field decreases the 
critical temperature because this introduces orbital frus¬ 
tration for the pairing fermions. (2) In the ERD case, 
at fixed non-zero Zeeman field, Tbkt increases relative 
to the case without spin-orbit coupling for all values of 
the binding energy. This is due to the emergence of a 
triplet component of the order parameter, induced by the 
presence of spin-orbit coupling. However, Tbkt never be¬ 
comes larger than the case of vanishing Zeeman and spin- 
orbit coupling fields, because of residual orbital effects. 
(3) The Clogston limit becomes infinite when spin-orbit 
coupling is present, in both the ERD case and the RO 
case. 

Finally, we constructed a 3D phase diagram, as a func- 
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T 


i-US-O 

US-2 

US-1 

i-US-O to US-2 
US-2 to US-1 
Tbkt 


FIG. 7. Three-dimensional phase diagram, as a function of the two-body binding energy E^, Zeeman held and temperature 
T, for v-r/vf = ud/Cf = 0.5. This phase diagram connects the zero-temperature phase diagram to the Berezinskii-Kosterlitz- 
Thouless critical temperature (red translucent surface). The zero-temperature phase diagram contains three uniform superhuid 
phases: i-US-O (red circles), US-2 (blue diamonds) and US-1 (green six-point stars). Finally, we have plotted the boundaries 
between the different US phases at non-zero temperature: i-US-O US-2 (blue surface) and US-2 —>■ US-1 (green surface). 


tion of the two-body binding energy, Zeeman field and 
temperature. We have extended the quantum phase 
transition lines to non-zero temperature, using the nodal 
structure of the quasiparticle excitation spectrum. The 
resulting phase diagram connects the zero-temperature 
result to the BKT-critical temperature, summarizing this 
paper. 
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Appendix A: Expanding the action np to qnadratic order in the phase 


In this section, we present the quadratic expansion of the effective action in more detail. As a first step, the second 
term in Eq. (1551) is expanded as 

Tr{ln[I + A-i (0,0)Fu(0, 90)]} « Tr[A-i (0, O)Fk(0, 80)] - ^TiiA-i (0,0)Fu(0, 90)A-i (0, O)Fk(0, 80)], (Al) 


leading to linear and quadratic terms in the expansion, which are treated separately. To calculate both these terms, 
we require the inverse of the matrix Ak.cj^ (0,0). Using symmetry relations, this inverse matrix can be written as a 
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function of only six elements 


C„(o>o) = 


D(k,uj,i) 


2 (k, -w„) 
AUO^) 

\ ^i,4(k,a;„) 


^l, 2 (k, (jJn) 
^2,2(k, UJn) 

“^l,4(ki 

A*Ak) 


^i,3(k) 

-AlAK^n) 

~^l,l(k, Wn) 

^ 1 . 2 (k, —LOn) 


Ai,4(k,a;„) \ 
^2,4(k) 
^l, 2 (k; Wn) 

2 (k,w„)y 


(A2) 


In this expression, the diagonal elements are equal to 

Ai,i(k, Wn) = {-iuJn + + hAi-iuJn - fk + hAi-iuJn - Ck - h^j - (-iWn + ^k + hA)\h±{\A? - |Ap(-ia;„ - ^k - hA), 
dl2.2 (k, U^n) — ( 'i^n A ^k hz) ( i^n ^k “f hz) ( iUJn ^k ) ( i^n “t" ^k hz )|Ii±(k)| |A| ( iuJn ^k “t" hz ) ■ 

(A3) 

Furthermore, the non-diagonal elements can be divided into those that depend both on the momentum k and the 
fermionic Matsubara frequency lOn 


Ai,2(k,a;„) = -h*j_{'k)\h±{'k)\'^ - |Ap/i;^(k) -b h*j_{'k){-iuJn - ■fk + hz){-iuJn - ■fk - hz), 
Ai,4(k,a;„) = |A|3 -b |hj_(k)p|A| - |A|(-za;„ -b fk + hz){-iu}n - Ck + hz), 

and those that only depend on the momentum k 

Ai,3(k) = -2/il(k)|A|(^k + /i.), 

A2,4(k) = 2hj.(k)|A|(ek-/i.). 

Finally, the determinant of the matrix Ak.(^„(0,0) is denoted by 

D{k,Uln) = (^-iujn + 4’^^(k)) -b 4”^!^)) “ 4"''^ (^)) “ 4~4k)) • 


(A4) 


(A5) 


(A6) 


Here, ep^^(k) denotes the quasiparticle energies. Through the rest of this derivation, we make use of the following 
symmetry properties of Il(k, a;„): 


£>(k,w„) = D*{k,u;n) = D(k,-uj„) = £)(-k,w„). 


(A7) 


We introduce the matrix Bk,i.,;„ {0, 09) = A^. ^ (0,0)Fk(0, 09) and then we write the linear term as 

Tr[Bk.c;„(6»,56l)] = ^ J dr J drJ2Yl D{k uj ) [Ai.i(k,a;n) + A 2 , 2 (k,a;n)] [ 4 (k) - £g_(k)] ^ , (A 8 ) 


where we used the definition £±(k) = — Cki with + |:[Vr(6>r,r)]^ and Ck = — Vr(0r,r) • k. This 

expression can be simplified by using boundary conditions in calculating the integrals over space and imaginary time. 
More specifically, we have 


4r^ = 0 . 

Ot 


(A9) 


Furthermore, in (I ASF we have used the fact that A*-(k,a;„) = Ai^i(k, — a;„) with i = {1,2}. By applying these 
properties, the linear expansion term can be simplihed to: 


Tr[Bk,..„(0,a0)] = 


2/3L2 


dr / J(k,O.„)[Vr(0r,r4 


k Uln 


where we have introduced 


A (k, iOfi 


D{k,UJn) 

Now, we take a look at the second contribution in the expansion m, which becomes 

1 


(AlO) 


(All) 


- iTr[Bk,.„(d,a0)Bk..„(d,ad)] =-^ JdTj 


k 


[£)(k,w„)]2 


(/(") (k, ojA + 7^'“^ (k, ojA + (k, . 

(A12) 
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(A13) 


We have divided the integrand in Eq. (IA12I) in three main terms, in order to keep track of this extensive expression. 
The terms in oj„) do not depend on spin-orbit coupling and are proportional to [e® (k)]^, [£^(k)]^ or (k)e® (k): 

/(^^(k,a;„) = [Ai i(k,a;„) + + 2Ai,2(k,u;„)A* 2(k,-w„)] {[£+(k)]2 [£?.(k)]^} 

+ 2 [|Ai.3(k)|2 + |A2,4(k)|2 + 2[Ai,4(k,a;„)]2] £t(k)£"_(k) 

Furthermore, the terms in /(^)(k,a;„) are proportional to [(^1)*]^ or with h^j_ = —— tct: 

7(")(k,c^„) = 2 [Al,{k,uj^) - Ai, 3 (k)A;, 4 (k)] [{hiy]\ 2{[Al,{k,u;nf - At,3(k)A2.4(k)} 

+ 4 [Ai,i(k,aj„)A2,2(k,w„) -|- |Ai_4(k,a;„)|^] 

Finally, the terms in w„) are proportional to the product £^(k)/i^ or £^ (k)(/.i)- 


/('''‘^(k,a;„) = 2 [£^(k)-7£?.(k)] 

X [Ai4(k,u;„)Ai,2(k,a;„) + Ai,3(k)Ai4(k,a;„) - Ai,2(k,a;„)A2,2(k, Wn) - Ai,4(k,a;„)A2_4(k)] 

+ 2 [£+(k) -I- £?.(k)] h^j_ 

X [-^1 2(k,-a;„)A2,2(k,w„) - Ai,4(k,u;„)A2,4(k) - Ai^i{k,u)n)Al 2(k, -w„) + A* 3(k)Ai,4(k, w„)] 

(A15) 

At this point, we retain only terms up to quadratic order in the phase field. Moreover, we use the fact that the 
Matsubara sum runs from —oo to oo to group terms together. Expression (IA13I) then becomes 


/^®)(k,a;„) = A(k,a;„) 


de^ 


dr 


+ S(k,a;„) 


86, 


dx 


dx dy 


86, 


8y 


kl 


(A16) 


where the following coefficients were defined: 


-A(k, OJn) — \ [—Af ^(k, LOn) + |Ai 3 (k)P — A| 2 (k, OOn) + |A 2 _ 4 (k)P — 2Ai 2 (k, Wji)A* 2 (k, —UJn) + 2A4^4(k, LUn)\ 

= 2 [A\ i(k,a;„) + |Ai,3(k)p -|- A\ 2(k, Wn) + |A2,4(k)p -I- 2Ay2(k, a;„)A 4 2(k, -uJn) + 2 A 4 4 (k,a;„)] 

(A17) 

Similarly, expression (IA14p can be written as 


/("•(k, „„) = C(k, % 5 (k. + £(k, „,J 


86,^r 

dy 


(A18) 


with the following coefficients: 

C(k, w„) = 4a^ [-A^_2(k,Wn) -f Ai, 3 (k)A 24 (k) -7 Ai,i(k,w„)A2,2(k, -f |Ai,4(k, u;„)p] 

P(k,a;n) = -4iaj [Af_2(k, Wn) - [A^ 2(k, w„)]^ - Ay3(k)A2_4(k) -|- A*_3(k)A2,4(k)] . (A 19 ) 

f (k, Wn) = 47^ [Af 2(k,a;„) - Ai,3(kjA^ 4(k) -|- Ai,i(k,w„)A2,2(k, w„) + |Ai7(k, u;„)p] 


In (|A19|) . we have made use of the fact that all terms proportional to k^ky vanish, because these are odd under a 
parity transformation of the integral over momentum k. Using the same reasoning, we can see that the momentum 
integral over 22(k,a;„) is equal to zero. Finally, we look at expression (IA15I) . which can be written as 

^ ^ (d6,.r\^ ^ ^d6,,rd0r.r , 77 ^, (d6,,r\^ 

/^®"Ak,u;„) = J^(k,a;„) ( 1 + 5(k,w„)-^-^-h 77(k, w„) ( 1 , (A20) 

with the following coefficients 

.F(k, w„) = -4iakxTa:(}^,^n) 

^(k,a;„) = 4[7A:,cry(k,a;„) - iakyT , (A21) 

'H(k,a;„) = 47A:yry (k, a;„) 

where ra;(k, u;„) is defined by 

r,;,(k,w„) =Ai,i(k,u;„)Ai, 2 (k,w„) - Ai,3(k)Ai7(k, a;„) - A* 2 (k,-a;„)A 2 , 2 (k, a;„) 

- Ai7(k,w„)A2,4(k) - Ai,i(k,u;„)A* 2 (k,-Wn) + A* 3 (k)Ai 7 (k,a;„) , 

+ Ai7(k, Wn)A2,2(k, UJn) + Ai^ 4 (k, a;n)A2 4(k) 
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and rj,(k, Wn) is defined by 


rj^(k,a;„) =Ai,i(k,a;„)Ai,2(k,w„) - Ai_3(k)Ai4(k, w„) + -Wri)^2,2(k, a;„) 

+ ^i,4(k, w„)A 24 (k) + Ai_i(k,a;„)^* 2(k, —uJn) — ^i,3(k)^i,4(k, w„) . 

+ ^l,2(k, W„)A2,2(k, iOn) + Ai_4(k, W„)A 2 4 (k) 


Using elementary algebra, we can show that the former two coefficients can be written as r 2 ,(k,a;„) oc kxf{k^,ky) 
and ry(k,a;„) oc fcy/(fc^, fc^), where / is a real function depending only on the square of the momentum components. 
Applying this form to expression (IA21I) . we see that expression t/(k, uJn) is proportional to k^ky and that its momentum 
integral is zero. Putting everything together, the total action can be written in a compact generic form: 





/^ p^y\ / d6/dx\ 

\dx dyj \Pyx Pvvj \de/dyj 


Here, we defined the coefficient 


A = 


1 

2/3L2 


^ A(k,w„), 

k,a;n 


which plays the role of the compressibility of the superfluid, and 

Pxx = X! (/3 + S(k,a;„)/c^ +C(k,a;„) +.F(k,w„) -^(k,a;„)) 

k,a;„ 

< Pxy — Pyx — A or 2 ^ ^ ^n'j^x^y ^n) “t“ ^(k, 

k,^;^ 

Pvv = X] (^ + ^+ ‘H{\^,UJn) - J (k, UJrS) 

k,a;n 


(A22) 


(A23) 


(A24) 


which play the role of the superfluid density tensor components pij. As mentioned previously, both the momentum 
integral over 22(k,a;„) and C/(k, w„) are zero. Furthermore, because all terms in H(k,a;„) are proportional to or 
ky, this part of the integral also vanishes due to the factor kxky. Hence, we have that pxy = Pyx = 0, making the 
superfluid density tensor pij diagonal. 

In the limit of spin-orbit coupling tending to zero, the action (jA22|l reduces to the known form (d^ . 















